# Code to support "Absolute and Relative Mobility: Two Frameworks for Connecting Intergenerational Mobility in Absolute and Relative Terms"
# (Sociological Methods and Research)

#####################
# NOTE TO USER:
####################

# Before running this file, ensure that the data file and function file are saved in user-set working directory

#####################
# Set up
####################

# Set working directory
# (User must set)
setwd("") 

# Read in data
dat <- read.csv("nlsy_dat.csv") 

# Install/load packages and set parallel computing parameters
list.of.packages <- c("Hmisc","dplyr","parallel","pbmcapply","mvtnorm","xtable")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)){
	install.packages(new.packages,dependencies=TRUE)
	}
lapply(list.of.packages, require, character.only = TRUE)
mcoptions <- list(preschedule=FALSE, set.seed=TRUE)
ncores <- detectCores()

# Source functions
source("functions.R")

####################
# Input parameters
####################

# Estimate parameters 
mu.p <- wtd.mean(dat$origin_inc_log,dat$weight)
mu.a <- wtd.mean(dat$dest_inc_log,dat$weight)
sigma.p <- sqrt(wtd.var(dat$origin_inc_log,dat$weight))
sigma.a <- sqrt(wtd.var(dat$dest_inc_log,dat$weight))
rho.reg <- coef(lm(I(scale(dest_inc_log))~I(scale(origin_inc_log)),weights=weight,data=dat))[2] 

# Create vectors of parameter values in the neighborhood of observed values 
adj.factor <- seq(from=.5,to=1.5,by=.05) 
mu.p.vec <- mu.p*adj.factor
mu.a.vec <- mu.a*adj.factor
sigma.p.vec <- sigma.p*adj.factor 
sigma.a.vec <- sigma.a*adj.factor
rho.reg.vec <- rho.reg*adj.factor 
rhoM.reg.vec <- (1-rho.reg)*adj.factor

####################
# Figure 1
####################

# Create input 
for(i in 1:(length(rhoM.reg.vec))) {
	j <- 5*(i-1)+50
	tmp <- a.pred.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec[i],dat$origin_inc_log,nrow(dat))
	out <- input.fig01.fun(tmp)
	assign(paste("input.fig01.bc",j,sep="."),out)
  }
input.fig01.a <- a.cond.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,0,dat$origin_inc_log) 

# Plot
pdf("fig01.pdf",height=4,width=12)
par(mfrow=c(1,3))
# (Panel A) 
  par(mar=c(5.5,5.5,4,2))
  plot(rhoM.reg.vec,
       apply(input.fig01.a$a.mob.mat[input.fig01.a$xp.above.mup==1 & input.fig01.a$xp.above.mua==1,],
             2,function(col) wtd.mean(x=col,weights=dat$weight[input.fig01.a$xp.above.mup==1 & input.fig01.a$xp.above.mua==1])),
       xlim=c(.2,.9),ylim=c(30,90),type="b",lty=1,pch=24,
       bty="n",xlab="1 - Intergenerational Correlation \n (Relative Mobility)",
       ylab="Mean % Upwardly Mobile \n (Absolute Mobility)",las=1,
       col=c(rep("black",10),"red",rep("black",10)),bg=c(rep("black",10),"red",rep("black",10)),
       main="% Experiencing Income Gain")
  points(rhoM.reg.vec,
         apply(input.fig01.a$a.mob.mat[input.fig01.a$xp.above.mup==1 & input.fig01.a$xp.above.mua==0,],
             2,function(col) wtd.mean(x=col,weights=dat$weight[input.fig01.a$xp.above.mup==1 & input.fig01.a$xp.above.mua==0])),
         type="b",lty=2,pch=18,
         col=c(rep("black",10),"red",rep("black",10)))    
  points(rhoM.reg.vec,
         apply(input.fig01.a$a.mob.mat[input.fig01.a$xp.above.mup==0 & input.fig01.a$xp.above.mua==0,],
             2,function(col) wtd.mean(x=col,weights=dat$weight[input.fig01.a$xp.above.mup==0 & input.fig01.a$xp.above.mua==0])),
         type="b",lty=3,pch=25,
         col=c(rep("black",10),"red",rep("black",10)),bg=c(rep("black",10),"red",rep("black",10)))         
  legend(.21,47,c("Par. Inc. Below Par. Mean","Par. Inc. Above Par. Mean \n and Below Adult Mean","Par. Inc. Above Both Means"),
         pch=c(25,18,24),col="black",pt.bg="black",cex=.7,bty="n")      
# (Panel B) 
  par(mar=c(5.5,5.5,4,2)) 
  plot(rhoM.reg.vec,
       c(input.fig01.bc.50[1,1],input.fig01.bc.55[1,1],input.fig01.bc.60[1,1],input.fig01.bc.65[1,1],
         input.fig01.bc.70[1,1],input.fig01.bc.75[1,1],input.fig01.bc.80[1,1],input.fig01.bc.85[1,1],
         input.fig01.bc.90[1,1],input.fig01.bc.95[1,1],input.fig01.bc.100[1,1],input.fig01.bc.105[1,1],
         input.fig01.bc.110[1,1],input.fig01.bc.115[1,1],input.fig01.bc.120[1,1],input.fig01.bc.125[1,1],
         input.fig01.bc.130[1,1],input.fig01.bc.135[1,1],input.fig01.bc.140[1,1],input.fig01.bc.145[1,1],
         input.fig01.bc.150[1,1]),
       xlim=c(.2,.9),ylim=c(-.5,1.5),type="b",lty=1,pch=16,
       bty="n",xlab="1 - Intergenerational Correlation \n (Relative Mobility)",
       ylab="Mean Expected Gain (Log Points) \n (Expected Absolute Mobility)",las=1,
       col=c(rep("black",10),"red",rep("black",10)),bg=c(rep("black",10),"red",rep("black",10)),
       main="Log Gain")
  points(rhoM.reg.vec,
       c(input.fig01.bc.50[1,2],input.fig01.bc.55[1,2],input.fig01.bc.60[1,2],input.fig01.bc.65[1,2],
         input.fig01.bc.70[1,2],input.fig01.bc.75[1,2],input.fig01.bc.80[1,2],input.fig01.bc.85[1,2],
         input.fig01.bc.90[1,2],input.fig01.bc.95[1,2],input.fig01.bc.100[1,2],input.fig01.bc.105[1,2],
         input.fig01.bc.110[1,2],input.fig01.bc.115[1,2],input.fig01.bc.120[1,2],input.fig01.bc.125[1,2],
         input.fig01.bc.130[1,2],input.fig01.bc.135[1,2],input.fig01.bc.140[1,2],input.fig01.bc.145[1,2],
         input.fig01.bc.150[1,2]),
         type="b",lty=2,pch=24,
         col=c(rep("black",10),"red",rep("black",10)),bg=c(rep("black",10),"red",rep("black",10)))   
  legend(.2,-.1,c("Below Mean Parental Inc.","Above Mean Parental Inc."),
         lty=c(1,2),pch=c(16,24),col="black",pt.bg="black",cex=.7,bty="n")       
# (Panel C) 
  par(mar=c(5.5,5.5,4,2)) 
  plot(rhoM.reg.vec,
       c(input.fig01.bc.50[2,1],input.fig01.bc.55[2,1],input.fig01.bc.60[2,1],input.fig01.bc.65[2,1],
         input.fig01.bc.70[2,1],input.fig01.bc.75[2,1],input.fig01.bc.80[2,1],input.fig01.bc.85[2,1],
         input.fig01.bc.90[2,1],input.fig01.bc.95[2,1],input.fig01.bc.100[2,1],input.fig01.bc.105[2,1],
         input.fig01.bc.110[2,1],input.fig01.bc.115[2,1],input.fig01.bc.120[2,1],input.fig01.bc.125[2,1],
         input.fig01.bc.130[2,1],input.fig01.bc.135[2,1],input.fig01.bc.140[2,1],input.fig01.bc.145[2,1],
         input.fig01.bc.150[2,1]),
       xlim=c(.2,.9),ylim=c(-.5,1.5),type="b",lty=1,pch=16,
       bty="n",xlab="1 - Intergenerational Correlation \n (Relative Mobility)",
       ylab="Mean Expected Gain (SD Units) \n (Expected Absolute Mobility)",las=1,
       col=c(rep("black",10),"red",rep("black",10)),bg=c(rep("black",10),"red",rep("black",10)),
       main="SD Gain")
  points(rhoM.reg.vec,
         c(input.fig01.bc.50[2,2],input.fig01.bc.55[2,2],input.fig01.bc.60[2,2],input.fig01.bc.65[2,2],
           input.fig01.bc.70[2,2],input.fig01.bc.75[2,2],input.fig01.bc.80[2,2],input.fig01.bc.85[2,2],
           input.fig01.bc.90[2,2],input.fig01.bc.95[2,2],input.fig01.bc.100[2,2],input.fig01.bc.105[2,2],
           input.fig01.bc.110[2,2],input.fig01.bc.115[2,2],input.fig01.bc.120[2,2],input.fig01.bc.125[2,2],
           input.fig01.bc.130[2,2],input.fig01.bc.135[2,2],input.fig01.bc.140[2,2],input.fig01.bc.145[2,2],
           input.fig01.bc.150[2,2]),
         type="b",lty=2,pch=24,
         col=c(rep("black",10),"red",rep("black",10)),bg=c(rep("black",10),"red",rep("black",10)))    
  legend(.2,-.1,c("Below Mean Parental Inc.","Above Mean Parental Inc."),
         lty=c(1,2),pch=c(16,24),col="black",pt.bg="black",cex=.7,bty="n")      
dev.off()

####################
# Figure 2
####################

# Create input
delta.vec <- c(-.5,-.25,0,.25,.5)
for(i in 1:(length(delta.vec))) {  
	out <- a.mua.pop.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,delta.vec[i],0)
	assign(paste("input.fig02",i,sep="."),out)
  }

# Plot
pdf("fig02.pdf",height=12)
par(mfrow=c(2,1))
# (Panel A)
  par(mar=c(5.5,5.5,4,2))
  plot(rhoM.reg.vec,input.fig02.1$a.mob.vec,
       xlim=c(.2,.9),ylim=c(20,100),type="b",lty=1,pch=19,
       bty="n",xlab="1 - Intergenerational Correlation \n (Relative Mobility)",
       ylab="% Upwardly Mobile \n (Absolute Mobility)",las=1,main="Upward Absolute Mobility",
       col=c(rep("black",10),"red",rep("black",10)))
  points(rhoM.reg.vec,input.fig02.2$a.mob.vec,
         type="b",lty=2,pch=18,
         col=c(rep("black",10),"red",rep("black",10)))    
  points(rhoM.reg.vec,input.fig02.3$a.mob.vec,
         type="b",lty=3,pch=8,cex=.7,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,input.fig02.4$a.mob.vec,
         type="b",lty=4,pch=20,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,input.fig02.5$a.mob.vec,
         type="b",lty=4,pch=17,cex=.7,
         col=c(rep("black",10),"red",rep("black",10)))   
  legend(.67,98,rev(c(expression(paste(delta," = -0.50")),expression(paste(delta," = -0.25")),
                      expression(paste(delta," = 0 (observed mean growth)")),
                      expression(paste(delta," = +0.25")),expression(paste(delta," = +0.50")))),
         lty=rev(c(1,2,3,4,4)),pch=rev(c(19,18,8,20,17)),col="black",cex=.7,bty="n")   
# (Panel B)
  plot(rhoM.reg.vec,input.fig02.1$mean.growth.pct.vec,
       xlim=c(.2,.9),ylim=c(-50,150),type="b",lty=1,pch=19,
       bty="n",xlab="1 - Intergenerational Correlation \n (Relative Mobility)",
       ylab="% Change in Mean Income \n (Mean Growth)",las=1,main="Mean Income Growth",
       col=c(rep("black",10),"red",rep("black",10)))
  points(rhoM.reg.vec,input.fig02.2$mean.growth.pct.vec,
         type="b",lty=2,pch=18,
         col=c(rep("black",10),"red",rep("black",10)))    
  points(rhoM.reg.vec,input.fig02.3$mean.growth.pct.vec,
         type="b",lty=3,pch=8,cex=.7,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,input.fig02.4$mean.growth.pct.vec,
         type="b",lty=4,pch=20,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,input.fig02.5$mean.growth.pct.vec,
         type="b",lty=4,pch=17,cex=.7,
         col=c(rep("black",10),"red",rep("black",10)))   
  legend(.67,138,rev(c(expression(paste(delta," = -0.50")),expression(paste(delta," = -0.25")),
                      expression(paste(delta," = 0 (observed mean growth)")),
                      expression(paste(delta," = +0.25")),expression(paste(delta," = +0.50")))),
         lty=rev(c(1,2,3,4,4)),pch=rev(c(19,18,8,20,17)),col="black",cex=.7,bty="n")   
dev.off() 

####################
# Figure 3
####################

pdf("fig03.pdf")
quads = function(colours=c("lightgray","darkgray","lightgray","darkgray")){
  limits = par()$usr
  rect(0.3,0.3,limits[2],limits[4],col=colours[1],border="NA")
  rect(-0.3,0.3,limits[1],limits[4],col=colours[2],border="NA")
  rect(-0.3,-0.3,limits[1],limits[3],col=colours[3],border="NA")
  rect(0.3,-0.3,limits[2],limits[3],col=colours[4],border="NA")
}  
x <- runif(20,-10,10) 
y <- runif(20,-10,10)
plot(x,y,type="n",bty="n",xlab="Relative Mobility",ylab="Absolute Mobility",
     xaxt="n",yaxt="n",xlim=c(-10,10),ylim=c(-10,10))
axis(1,labels=0,at=0,lwd=2,lty=1)
axis(2,labels=0,at=0,lwd=2,lty=1)
abline(h=0,lwd=2)
abline(v=0,lwd=2)
quads()     
text(5+.3,5,"Upwardly \n Mobile",font=2)  
text(-5-.3,-5,"Downwardly \n Mobile",font=2) 
text(-5-.3,5,"Floating with \n Rising Tide",font=2)
text(5+.3,-5,"Emerging from \n Retreating Tide",font=2)  
dev.off()

####################
# Figure 4
####################

# Create input 
fig04.point.above <- fig04.fun.all(dat[dat$origin_inc_above_mean==1,])
fig04.point.below <- fig04.fun.all(dat[dat$origin_inc_above_mean==0,])
nboots <- 1000  
fig04.boot.out.above <- fig04.fun.boot.pl(dat,1,nboots)
fig04.boot.out.below <- fig04.fun.boot.pl(dat,0,nboots)
fig04.cil.above <- fig04.fun.boot.collect(fig04.boot.out.above,nboots)$j.tab.1.q025 
fig04.ciu.above <- fig04.fun.boot.collect(fig04.boot.out.above,nboots)$j.tab.1.q975 
fig04.cil.below <- fig04.fun.boot.collect(fig04.boot.out.below,nboots)$j.tab.1.q025 
fig04.ciu.below <- fig04.fun.boot.collect(fig04.boot.out.below,nboots)$j.tab.1.q975 

# Plot 
pdf("fig04.pdf",width=14) 
mat.plot <- matrix(c(1,2,2),nrow=1,ncol=3,byrow = TRUE)
layout(mat=mat.plot)
# (Panel A) 
options(scipen=5)
par(mar=c(5.5,5.5,4,3)) 
plot(100*(dat$dest_rank-dat$origin_rank),dat$dest_inc_dol-dat$origin_inc_dol,
     bty="n",type="p",pch=16,col=gray(.5,.25),ylim=c(-170000,330000),
     las=1,ylab="Absolute Mobility ($) \n ",xlab="Relative Mobility (Ranks)",
     cex.lab=1.2,
     main="Absolute and Relative Mobility:\nJoint Experiences")
abline(h=0,lty=3,col="darkgray")
abline(v=0,lty=3,col="darkgray") 
text(10*5.2,330000,bquote(paste("Upwardly Mobile: ",
     .(round(100*wtd.mean(dat$dest_rank-dat$origin_rank>0 & dat$dest_inc_dol-dat$origin_inc_dol>0,dat$weight,na.rm=T),0)),"%"))) 
text(10*-5.2,-160000,bquote(paste("Downwardly Mobile: ",
     .(round(100*wtd.mean(dat$dest_rank-dat$origin_rank<0 & dat$dest_inc_dol-dat$origin_inc_dol<0,dat$weight,na.rm=T),0)),"%")))
text(10*-5.2,330000,bquote(paste("Floating with Rising Tide: ",
     .(round(100*wtd.mean(dat$dest_rank-dat$origin_rank<0 & dat$dest_inc_dol-dat$origin_inc_dol>0,dat$weight,na.rm=T),0)),"%")))
text(10*5.2,-160000,bquote(paste("Emerging from Retreating Tide: ",
     .(round(100*wtd.mean(dat$dest_rank-dat$origin_rank>0 & dat$dest_inc_dol-dat$origin_inc_dol<0,dat$weight,na.rm=T),0)),"%"))) 
# (Panel B)
var.names <- rev(c("White - Black","","","White - Hispanic","","","White - Other","","","Black - Hispanic","","","Black - Other","","",
                   "Stable Two Par. - Not","","","Native-Born Par. - Not","","")) 
cols.fig04 <- rev(rep(c("black","black","black"),length(var.names)))
pch.fig04 <- rev(rep(c(24,25,21),length(var.names)/3))
lty.fig04 <- rev(rep(c(1,1,2),length(var.names)/3))   
  par(mar = c(5,8.3,4,2))
  plot(NULL,                              
       xlim = c(-40, 130),                      
       ylim = c(0, length(var.names)/2+3), 
       axes = F, xlab = NA, ylab = NA, bty = "n",
       main="Absolute and Relative Mobility:\nGroup Differences in Joint Experiences")   
  j <- .1
  j.at <- .1
  for (i in 1:length(var.names)) {      
    points(rev(t(fig04.point.below))[i], j, cex = .5, col = cols.fig04[i], bg = cols.fig04[i], pch = pch.fig04[i])
    lines(c(rev(t(fig04.ciu.below))[i], rev(t(fig04.cil.below))[i]), c(j, j), col = cols.fig04[i], lty = lty.fig04[i])  
    points(90+rev(t(fig04.point.above))[i], j, cex = .5, col = cols.fig04[i], bg = cols.fig04[i], pch = pch.fig04[i])
    lines(c(90+rev(t(fig04.ciu.above))[i], 90+rev(t(fig04.cil.above))[i]), c(j, j), col = cols.fig04[i], lty = lty.fig04[i]) 
    j <- ifelse(i==3 | i==6 | i==9 | i==12 | i==15 | i==18 | i==21,j+1,j+.4) 
    j.at <- c(j.at,j)    
  } 
  abline(h=mean(j.at[3:4]),col="darkgray",lty=3)
  abline(h=mean(j.at[6:7]),col="darkgray",lty=3)
  abline(h=mean(j.at[9:10]),col="darkgray",lty=3)
  abline(h=mean(j.at[12:13]),col="darkgray",lty=3)
  abline(h=mean(j.at[15:16]),col="darkgray",lty=3)
  abline(h=mean(j.at[18:19]),col="darkgray",lty=3)
  abline(h=mean(j.at[21:22]),col="darkgray",lty=3)
  abline(v=0,col="darkgray",lty=2)
  abline(v=50,col="darkgray",lty=1)
  abline(v=90,col="darkgray",lty=2)
  axis(side = 1, labels = c(-40,0,40,-40,0,40), at = c(-40,0,40,50,90,130),cex.axis=1)
  axis(side = 2, labels = var.names, at = j.at[-length(j.at)],las=1,cex.axis=1)
  legend(109,13.2,c("Upward","Downward","Floating"),bty="n",cex=1,pt.cex=.5,
         lty=rev(lty.fig04),pch=rev(pch.fig04),pt.bg=rev(cols.fig04))
  mtext(c("Section A","Section B","Section C"),
          side = 2, line = 4, cex = .7, las = 1, at = j.at[c(21,6,3)]+.8,font=2)
  mtext("Percentage Point Difference", side = 1, line = 3, cex = .8, at = 0)
  mtext("Percentage Point Difference", side = 1, line = 3, cex = .8, at = 90)
  text(0-.1,13.5,"Below Origin Mean Income",font=2,cex=1)
  text(90-.1,13.5,"Above Origin Mean Income",font=2,cex=1)
dev.off()

####################
# Figure A1
####################

# Create input
n <- 1e6
t.sim.d3 <- t.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,3)
  t.sim.d5 <- t.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,5)
  t.sim.d10 <- t.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,10)
  t.sim.d30 <- t.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,30)
g.sim.s1 <- gamma.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,1)
  g.sim.s5 <- gamma.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,5)
  g.sim.s10 <- gamma.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,10)
  g.sim.s100 <- gamma.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,100)
  g.sim.s1000 <- gamma.sim.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,n,1000)
a.from.r <- a.pop.fun(mu.a,mu.p,sigma.a,sigma.p,1-rhoM.reg.vec,0) 

# Plot
pdf("fig_a01.pdf", width=14)
par(mfrow=c(1,2))
# (Panel A)
par(mar=c(5.5,5.5,4,2))
plot(rhoM.reg.vec,t.sim.d3$pct.up,
       xlim=c(.2,.9),ylim=c(30,90),type="b",lty=1,pch=19,
       bty="n",xlab="1 - Intergenerational Correlation \n (Relative Mobility)",
       ylab="% Upwardly Mobile \n (Absolute Mobility)",las=1,
       col=c(rep("black",10),"red",rep("black",10)),
       main="Assumed Distribution: \nBivariate t")
  points(rhoM.reg.vec,t.sim.d5$pct.up,
         type="b",lty=2,pch=18,
         col=c(rep("black",10),"red",rep("black",10)))    
  points(rhoM.reg.vec,t.sim.d10$pct.up,
         type="b",lty=3,pch=17,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,t.sim.d30$pct.up,
         type="b",lty=4,pch=20,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,a.from.r$a.mob,
         type="b",lty=4,pch=8,cex=.6,
         col=c(rep("black",10),"red",rep("black",10)))   
  legend(.75,90,c("DF = 3","DF = 5","DF = 10","DF = 30","Normal"),
         lty=c(1,2,3,4,4),pch=c(19,18,17,20,8),col="black",cex=.7,bty="n")      
# (Panel B) 
par(mar=c(5.5,5.5,4,2))
plot(rhoM.reg.vec,g.sim.s1000$pct.up,
       xlim=c(.2,.9),ylim=c(30,90),type="b",lty=1,pch=19,
       bty="n",xlab="1 - Intergenerational Correlation \n (Relative Mobility)",
       ylab="% Upwardly Mobile \n (Absolute Mobility)",las=1,
       col=c(rep("black",10),"red",rep("black",10)),
       main="Assumed Distribution: \nBivariate Gamma")
  points(rhoM.reg.vec,g.sim.s100$pct.up,
         type="b",lty=2,pch=18,
         col=c(rep("black",10),"red",rep("black",10)))    
  points(rhoM.reg.vec,g.sim.s10$pct.up,
         type="b",lty=3,pch=17,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,g.sim.s5$pct.up,
         type="b",lty=3,pch=15,
         col=c(rep("black",10),"red",rep("black",10))) 
  points(rhoM.reg.vec,g.sim.s1$pct.up,
         type="b",lty=4,pch=20,
         col=c(rep("black",10),"red",rep("black",10)))  
  points(rhoM.reg.vec,a.from.r$a.mob,
         type="b",lty=4,pch=8,cex=.6,
         col=c(rep("black",10),"red",rep("black",10)))   
  legend(.72,90,rev(c("Shape = Obs./1000","Shape = Obs./100","Shape = Obs./10","Shape = Obs./5","Shape = Obs.","Normal")),
         lty=rev(c(1,2,3,3,4,4)),pch=rev(c(19,18,17,15,20,8)),col="black",cex=.7,bty="n")     
dev.off() 

####################
# Table 1
####################

a.pred <- a.pred.fun(mu.a,mu.p,sigma.a,sigma.p,rho.reg,dat$origin_inc_log,nrow(dat))
tab01 <- round(tab01.fun(a.pred,dat$origin_inc_above_mean,dat$weight),2)
nboots <- 1000  
tab01.boot.out <- tab01.fun.boot.pl(dat,nboots)
tab01.cis <- tab01.fun.boot.collect(tab01.boot.out,nboots)$j.tab.1.ci
tab01.with.cis <- rbind(tab01[1,],tab01.cis[1,],tab01[2,],tab01.cis[2,],tab01[3,],tab01.cis[3,],tab01[4,],tab01.cis[4,])
rownames(tab01.with.cis) <- c(rownames(tab01)[1],"CI1",rownames(tab01)[2],"CI2",rownames(tab01)[3],"CI3",rownames(tab01)[4],"CI4")
xtable(tab01.with.cis)

####################
# Table 2
####################

tab02 <- round(tab02.fun(dat),2)
nboots <- 1000   
tab02.boot.out <- tab02.fun.boot.pl(dat,nboots)
tab02.boot.coll <- tab02.fun.boot.collect(tab02.boot.out,nboots) 
tab02.with.cis <- tab02.fun.boot.combine(tab02,tab02.boot.coll,1)
xtable(tab02.with.cis) 

####################
# Table 3
####################

tab03 <- round(tab03.fun.all(dat),2)
nboots <- 1000   
tab03.boot.out <- tab03.fun.boot.pl(dat,nboots)
tab03.cis <- tab03.fun.boot.collect(tab03.boot.out,nboots)$j.tab.1.ci 
tab03.with.cis <- rbind(tab03[1,],tab03.cis[1,],tab03[2,],tab03.cis[2,],tab03[3,],tab03.cis[3,],tab03[4,],tab03.cis[4,],
                        tab03[5,],tab03.cis[5,],tab03[6,],tab03.cis[6,],tab03[7,],tab03.cis[7,])
rownames(tab03.with.cis) <- c(rownames(tab03)[1],"CI1",rownames(tab03)[2],"CI2",rownames(tab03)[3],"CI3",rownames(tab03)[4],"CI4",
                              rownames(tab03)[5],"CI5",rownames(tab03)[6],"CI6",rownames(tab03)[7],"CI7")
xtable(tab03.with.cis)

####################
# Table A1
####################

tab.a01 <- rbind(c(mu.a,sigma.a),c(mu.p,sigma.p), 
                 c(100*wtd.mean(dat$white,dat$weight,na.rm=T),999),
                 c(100*wtd.mean(dat$black,dat$weight,na.rm=T),999),
                 c(100*wtd.mean(dat$hispanic,dat$weight,na.rm=T),999),
                 c(100*wtd.mean(dat$other,dat$weight,na.rm=T),999),
                 c(100*wtd.mean(dat$stable_par,dat$weight,na.rm=T),999),
                 c(100*wtd.mean(dat$native_par,dat$weight,na.rm=T),999),
                  c(nrow(dat),999))
colnames(tab.a01) <- c("Mean","SD")
rownames(tab.a01) <- c("Adult income","Parental income","% White","% Black","% Hispanic","% Other",
                        "% Stable two-parent childhood family structure","% Native-born parents","N")
xtable(tab.a01,digits=2)
